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Abstract. 

The equation of state (EOS) is an important issue in numerical simulation codes 
for relativistic hydrodynamics. We describe a code for adiabatic flows, employing an 
EOS which is simple and yet approximates very closely the EOS of perfect gas in 
relativistic regime. We also describe a code for isothermal flows, where the EoS is 
trivially given. 



1. Introduction 

Highly energetic phenomena, which are relativistic in nature, are c ommon in as tro- 



physical environments: accretion disks around black holes (see, e.g., [M iller 2007, for 



review), relativistic jets from Galactic sources (see, e.g.. iMirabel & Rodrig uez., 19991 



for review), extragalactic jets from active galac tic nuclei ( see, e.g., Zensus 1 19971. for 



review), and gamma-ray bursts (see, e.g., iMeszarosI l2002l for review). Gas in such 



relativistic phenomena is characterized by its relativistic fluid speed (v — > c) and/or 
relativistic sound speed (c^ — > c/ V3). 

Numerical codes for relativistic hydrodynamics (RHDs) have been successfully 
built, based on schemes that were originally developed for codes for non-relativistic 
hydrodynamics. Most codes employed the equation of state (EoS) of the gas, which 
was designed for the gas with a constant ratio of specific heats and so is essentially 
valid only for the ga s in either n on-relativistic or ultra-relativistic regime (see, e.g.. 



Martf & Mulled 120031 : IWilson & Mathews .2003,, for reviews). The correct EoS (see 



Section 3), however, involves the specific enthalpy expressed in terms of the modified 



Bessel function s (seelSy nge 1957). While codes employing th e correct EoS has been 
built (see, e.g.. iFalle & Komissarov 19961 : IScheck et al.ll2002l) . they normally comes 



with an extra cost of computation time. On other othe r hand, approxi r tiate EoSs that 
mimic the corre ct EoS have been suggested (see, e.g.. lMathewslll97fl : IServic 3 119861 : 



Ryu et al.l 120061) . an d recently, code s emplo ying tho se approximate EoSs have been 



introduced (see, e.g.. lMignone et al.i r2005: Ry u et al.l l2006) 



In this paper, we describe a RHD code for adiabatic flows, which was presented 
in lRvu et all (120061) : an EoS, which is simple and an algebraic function of temperature, 
was employed. We also describe a new RHD code for isothermal flows; the meaning 
of isothermality in relativistic regime is discussed. The steps necessary to build codes 
including the transformation from the conserved quantities to the primitive quantities 
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and the eigen-structure are presented. Finally, shock tube tests performed with codes 
based on the Total Variation Diminishing (TVD) scheme are presented. 



2. Relativistic Hydrodynamic Equations 

The special RHD equations for an ideal fluid in the laboratory frame of reference can 
be written as 

dD 5 / X 

dMi d I X 

+ ^(^<v, + p<5,;)-0, (2) 



dt dxj 
dE_ d_ 



f ,^[(£,^),.] = 0, (3) 



where D, Mi, and E are th e mass density, momentu m density, and total energy den- 
sity, respectively (see, e.g. jLandau & Lifshitzlll959l) . The conserved quantities in the 
laboratory frame are expressed as 

D = Fp, Mi = Y^phvi, E = T^ph - p, (4) 

where p, v,-, p, and h are the proper mass density, fluid three-velocity, isotropic gas 
pressure and specific enthalpy, respectively, and the Lorentz factor is given by 

F = ^ with = vl + vl + v\. (5) 

Here, h = {e + p)lp, where e is the sum of the proper internal and rest-mass energy 
densities. In the above, the Latin indices (e.g., /) represents spatial coordinates and the 
conventional Einstein summation is used. The speed of light is set to unity (c = 1) 
throughout this paper. 



3. Code for Adiabatic Flows 

Equations (1) - (3) with the EoS, h = h{p,p), form a hyperbolic set of conservation 
equations for adiabatic flows. The correct EoS for th e single-corn ponent perfect gas in 
relativistic regime (hereafter RP) can be derived (see ISynge|[l957l) : it is given as 

= (6) 

where K2 and A'3 are the modified Bessel functions of the second kind of order two and 
three, respectively. Here, & = p/pis a. temperature-like variable. Using the EoS in (6), 
however, poses a difficulty, because the inverse, that is, as a function of h, can not be 
expr essed as a s imple form. 

iRvu et aP (l2006l) introduced an approximate EoS (hereafter RC), 

30-^2 ^ ^ 
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Figure 1 . Polytropic index, n, and sound speed, Cj, as a function of temperature, 
= p/p, for adiabatic flows employing different EoSs. 



RC mimics very closely RP. Figure 1 compares the polytropic index and sound speed 

(8) 



dh 

n=P^ - 1, 
op 



2 p dh 
* nhdp' 



for RP and RC as well as for the EoS with a constant ratio of specific heats, y, (hereafter 
ID, standing for ideal gas). The specific enthalpy h for RC fits h for RP within the error 
of 0.8 %. 

Building codes based on upwind schemes requires the eigen-structure (eigenvalues 
and eigenvectors) for the relevant hyperbolic set of conservation equations. The eigen 



structu re for RHD equations for adiabatic flows in (1) - (3) are given in IRvu et al. 
(I2006h . The eigenvalues are 



a5 



(l-c2)v,-c,/r- VG 

a2 = ai, - a4 = Vx, 

{i-cl)vx + cjr- VG 

1 - c^y ' 



(9) 



(10) 



(11) 



where 2 = 1 - - c^(Vy + v^), assu ming that the flow varies along the A-direction. For 
the right and left eigenvectors, refer IRvu et al.l(l2006h . 

Equations (1) - (3) evolve the conserved quantities D, Mj, and E, but the primitive 
quantities, p, v,, and p, are necessary to calculate the eigenvalues and eigenvectors. 




Combining (4) - (5) along with (7) results in 

MVr2 - 1 [3£r(8r2 - i) + 2D(i - Ar^)] 

= [4(M2 + E^)T^ - (M2 + 4£2)] _ 2D(4£r - D){Y^ - 1). (12) 

Further simplification reduces this to an equation involving the 8* power of Y. The 
above can be solved numerically to get v. And then, the rest of the primitive quantities 
can be solved. 

A RHD code for adiabatic flows was built based on the TVD scheme (IRvu et all 



20061) . Figure 2 show a shock tube test comparing the results with RC and ID: initially 
PL ^ PR ^ 1, PL = 10^ PR ^ 10"2, Vp,L = Vp^R = 0, v,,L = 0.9, v,^R = 0.99, and 
^end - 0.75. Here, the subscripts L and R denote the quantities in the left and right 
states of the initial discontinuity at x = 0.5, and fend is the time when the solutions 
are presented. And Vp and Vt are the velocity components parallel and transverse to 
the propagation of structures (i.e., the x-direction). The ID solution with y = 4/3 
matches well to the RC solution in the left of the contact discontinuity where the flow 
has » 1. But a difference is obvious in the region between contact discontinuity and 
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Figure 3. A relativistic shock tube from the code for isothermal flows and from 
the code for adiabatic flows employing RC with different initial densities, p - I, 1/3, 
1/10. 



shock because ~ 1 there. The ID solution is clearly different from the RC solution, 
indicating the importance of using correct EoSs. 



4. Code for Isothermal Flows 



As in non-relativistic hydrodynamics (see, e.g.. lKim et alJll999h . a code for isothermal 
flows, where the EoS is given by p = e with a constant sound speed c^, can be 
built in RHDs. Such EoS arises in several important situations. Most of all, when the 
constituent particles are ultra-relativistic or the fluid is dominated by radiation, that is, 
e » p, the sound speed goes Cs — > 1/ V3 and the EoS becomes p = (l/3)e. Also the 
EoS in degenerate m atters, such as in white dwarfs and neutron stars, may be modeled 
as p = c^^ e (see, e.g. IWeinbergll 19721) . 

Equations (2) - (3) with p = c^^ e form a complete, hyperbolic set of conservation 
equations for isothermal flows. The eigenvalues have the same form as those for adia- 
batic flows in (9) - (1 1), except two instead of three; this is expect by considering the 
nature of weak solutions of initial value problems in RHDs (see, e.g. iSmoUer & Temple 
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19931) . The left and right eigenvectors are substantially simpler than those for adiabatic 



flows. (They will be published elsewhere due to the page limit of this proceeding pa- 
per.) The calculation of the primitive quantities from the conserved quantities is also 
substantially simpler; 



£-2 - 4c?M2 + (1 + c2)£ M 

v = ^^ — , p = — -E. (13) 

2c?M V 

Each component of velocity can be calculated with v, = (Mi/M)v. Again a RHD code 
for isothermal flows was built based on the TVD scheme. Figure 3 shows a shock tube 
test comparing the results from the code for isothermal flows and the code for adiabatic 
flows with RC: initially pi - I, pR - 0.1, Vp^ - 0.1, v^r = 0, - 0, v,^r = 0.9, 
and fend = 0.75. For the code for adiabatic flows, pt = Pr = 1, 1/3, and 1/10. The 
adiabatic solution approaches the isothermal solution, as p ^ small or ^ large. Our 
test indicates that for > a few xlO, the two solutions for isothermal and adiabatic 
flows become indistinguishable. 

In principle, isothermal flows can be simulated with codes for adiabatic flows, as 
well. However, there are advantages of using codes for isothermal flows: 1) Codes for 
isothermal flows are faster than those for adiabatic flows, because the eigenvectors are 
simpler. With our codes based on the TVD scheme, the isothermal version is about 
1.5 to 2 times faster than the adiabatic version. 2) Codes for isothermal flows should 
be numerically more robust than those for adiabatic flows, because one mode, which 
is the entropy mode, is less. Our tests have shown that shocks and discontinuities of 
Vf are better resolved with the code for isothermal flows. In addition, we expect that 
the numerical dissipation would be smaller in the code for isothermal flows (not shown 
here). 
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